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ABSTRACT 

We simulate fragmentation and gravitational collapse of cold, magnetized molecular clouds. 
We explore the nonlinear development of an instability mediated by ambipolar diffusion, in 
which the collapse rate is intermediate to fast gravitational collapse and slow quasistatic collapse. 
Initially uniform stable clouds fragment into elongated clumps with masses largely determined 
by the cloud temperature, but substantially larger than the thermal Jeans mass. The clumps are 
asymmetric, with significant rotation and vorticity, and lose magnetic flux as they collapse. The 
clump shapes, intermediate collapse rates, and infall profiles may help explain observations not 
easily fit by contemporary slow or rapid collapse models. 
Introduction 



The interstellar magnetic field plays an impor- 
tant role in the dynamics of molecular clouds and 
the collapse of dense cloud cores into protostars. 
Magnetic pressure and tension combine with ther- 
mal and turbulent kinetic pressure to resist gravi- 
tational collapse. The role of the magnetic field 
is often simply characterized by a critical mass 
Merit which depends on the magnetic flux thread- 
ing the cloud (Mestel & Spitzer 1956; Mouschovias 
& Spitzer 1976; Tomisaka, Ikeuchi, & Nakamura 
1988; McKee et al. 1993; Mestel 1999). Clouds 
with M > Merit are termed supercritical and col- 
lapse on a dynamical time-scale: the magnetic 
field can have a moderate effect on the morphol- 
ogy of collapse and can slow collapse in the cloud 
envelope {e.g. Black & Scott 1982), but can- 
not significantly slow collapse in the core. Clouds 
with M < Merit are termed subcritical, and evolve 
on a longer time-scale as magnetic support is lost 
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due to ambipolar diffusion. The magnetic field 
is redistributed within the cloud so that the in- 
ner parts become supercritical. The cloud is then 
differentiated into a dynamically collapsing core 
with a magnetically supported envelope (Ciolek & 
Mouschovias 1993; Fiedler & Mouschovias 1993; 
Basu & Mouschovias 1994; Ciolek & Mouschovias 
1994; Safier, McKee, & Stabler 1997; Ciolek & 
Konigl 1998). Much progress has been made in 
following this type of evolution through 6 or more 
orders of magnitude of increase in central den- 
sity, including the effects of rotation as well as de- 
tailed chemistry and grain physics. The outcomes 
of these calculations include detailed density, bulk 
velocity, ion-neutral drift velocity, magnetic field, 
and grain and ion abundance profiles in axisym- 
metric clouds, as well as an appreciation of the 
timescales, rate of magnetic flux loss, and role of 
magnetic braking in this mode of isolated star for- 
mation. 

This theoretical picture can be observationally 
tested (see recent reviews by Evans 1999; Myers, 
Evans, & Ohashi 1999). So far, the results are am- 
biguous. Flow toward an isolated infrared source 
in the Bok globule B335 is well described by an 
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inside-out collapse model (Zhou et al. 1990, 1993; 
Zhou 1995). The density distribution and mea- 
sured magnetic fieldstrength in the cloud Bl have 
been fit by a model with a subcritical envelope and 
a core which has evolved to a supercritical state 
by ambipolar drift (Crutcher et al. 1994). On the 
other hand, in some respects the existing models 
appear to be incomplete. Clumps and cores are 
not axisymmetric; Myers et al. (1991) surveyed 16 
dense cores a few tenths of a parscc in size and 
pointed out that at least 6 of them are likely to 
be prolate. Ryden (1996) made a statistical argu- 
ment, based on a larger sample, that clumps and 
globules are more likely prolate or triaxial than 
oblate. Ward-Thompson, Motte, & Andre (1999) 
showed that asphericity appears also at smaller 
scales. Collapse guided by a magnetic field could 
produce oblate clouds, but not prolate ones. It 
also appears unlikely that rotation accounts for 
the flattening; this has been shown quantitatively 
in the case of L1527 (Ohashi et al. 1997). Thus, 
the shapes are unexplained. Moreover, in some 
cases in which infall has been measured directly, 
it is more spatially extended, with faster veloc- 
ities in the outer parts, than expected from the 
standard models of inside-out collapse or gravi- 
tational motion driven by ambipolar drift. This 
has been shown in the case of LI 544 by Tafalla 
et al. (1998), and for 6 other starless cloud cores 
by Gregersen (1998). (Interestingly, a model of 
LI 544 has been recently constructed by Ciolek & 
Basu (2000) to match the observations of Tafalla et 
al. (1998). The model requires a mass-to-flux ra- 
tio more nearly critical than previously published 
models by the same authors, and appears to dis- 
play the same intermediate collapse discussed in 
this paper, although the authors do not call it out 
as such.) Finally, there are observations which 
relate to the timescale for collapse of molecular 
clouds into protostars. Lee & Myers (1999) find 
collapse timescales of ~ 0.3-1.6 Myr from the ra- 
tio of the numbers of starless cores to cores with 
embedded young stellar objects. They state that 
this requires collapse 2-44 times faster than am- 
bipolar drift models. These observations suggest 
that another ingredient may be required to explain 
the collapse of molecular clouds: the decay of tur- 
bulent support (Myers & Lazarian 1998), or, as 
we explore in this paper, a magnetogravitational 
instability mediated by ambipolar drift. 



The magnetic field likely also plays an impor- 
tant role in other aspects of collapse, which are 
still incompletely understood. The mechanism by 
which molecular clouds fragment, and the masses 
and morphology of those fragments, is clearly of 
importance to the stellar initial mass function and 
the origin of binary systems. The magnetic field 
can exert strong forces on many scales, affecting 
fragmentation {e.g. Boss 1997, 1999). The fleld 
may also be a source of kinetic energy in the cloud, 
if the free energy of an ordered field can be released 
as turbulent motions. These issues have not been 
studied thoroughly in the detailed axisymmetric 
models of gravitational contraction with ambipo- 
lar drift referenced above because they would re- 
quire calculations with no constraints on spatial 
symmetry. 

In this paper, we address some of these issues by 
studying a simple problem: the evolution of small 
perturbations to an initially uniform, magnetically 
subcritical sheet of weakly ionized gas with a uni- 
form magnetic field perpendicular to its plane. 
The perturbations evolve under the influence of 
magnetic tension, self gravity, thermal pressure, 
and ambipolar drift. Typically the sheet breaks 
up into a small number of fragments of elongated 
shape which are collapsing, losing magnetic flux 
through ambipolar drift, and interacting gravita- 
tionally with one another. The magnetic fields as- 
sociated with these asymmetric clumps generate 
local vorticity (the net angular momentum of the 
sheet is identically zero). The characteristic vor- 
ticity structure is a vortex pair which flanks each 
clump and is associated with strong streaming mo- 
tions along it. The clump masses are typically of 
order 1-10 M©, and scale with temperature like 
the Jeans mass, but are typically larger because; 
of magnetic siipport. The main features of col- 
lapse in this geometry were predicted by the linear 
stability analysis of Zweibel (1998; hereafter Z98); 
there is also some overlap with the earlier stability 
analysis in 3D by Langer (1978). Here we follow 
the evolution into the nonlinc;ar rc;ginie and follow 
the growth of density fluctuations from .01 to up 
to 10 times the mean surface density. Collapse oc- 
curs on an intermediate time-scale, slower than 
the dynamical or free-fall time-scale, but faster 
than the ambipolar diffusion time-scale. As noted 
above, this intermediate collapse rate may be ob- 
served in some clouds {e.g. Lee & Myers 1999; 
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Grcgcrscn 1998) . As no restrictions arc placed on 
the morphology of the clumps, we can begin to ex- 
plore the nature of flux loss and collapse in more 
complicated geometries than isolated axisymmet- 
ric clouds. One outcome suggested by linear the- 
ory which we have not resolved is whether stored 
magnetic energy is converted to turbulence, as 
there is no stored magnetic energy in the system 
initially. 

The unperturbed initial geometry, governing 
equations, and linear theory are described in §2. 
Section 3 contains a description of the numerical 
method and main results: the collapse rate is dis- 
cussed in §3.1, the relationship between magnetic 
field B and density p or surface density cr in §3.2, 
size of fragments in §3.3, velocity structure in §3.4, 
and distribution of energy in §3.5. The validity 
of the approximations is discussed in §4, and the 
summary and conclusions are in §5. 

2. Governing Equations and Linecir The- 
ory 

We simulate a flat slab of cold gas with the mag- 
netic field initially perpendicular to the slab (the 
z direction). Z98 discusses the linear theory for 
this model when the cloud temperature T equals 
0. In this section we review the model in the more 
general case T>0. We recall the results of the lin- 
ear theory and discuss the consequences of adding 
thermal pressure. 

A flat slab-like model has observational and 
theoretical motivation: molecular clouds com- 
monly have sheet or filamcnt-like structure (al- 
though detailed, high-resolution information on 
the fleld orientation in such structures is not yet 
available for many objects). In a fairly quiescent 
environment, a roughly spherical molecular cloud 
with a large-scale, dynamically significant, ordered 
magnetic field will relax into a pancake or slab 
as matter drains down the field lines. Magnetic 
forces will allow comparatively little contraction 
perpendicular to the field direction, resulting in 
a slab with a predominantly perpendicular field. 
Such slabs could also be formed by shock waves 
propagating parallel to the local magnetic field. 

We use a simple initial state with small (<1%) 
perturbations. The boundary conditions are pe- 
riodic in the horizontal (x and y) directions and 
the initial surface density (tq and magnetic field 



Bzo are uniform. The unperturbed initial state 
was chosen for computational simplicity, as self- 
consistent finite disk and slab-like equilibrium 
states cannot generally be described by simple 
expressions in closed form. {e.g. Parker 1974; 
Mouschovias 1976; Baureis et al. 1989; Mestel & 
Ray 1985). We note that our initial state is not 
technically one of static equilibrium, but rather 
a version of the commonly used Jeans swindle, 
in which the unperturbed gravitational poten- 
tial is discarded (see e.g. discussion in Binney 
& Tremaine (1987)). 

2.1. Governing Equations 

We begin with the equations of ideal magne- 
tohydrodynamics for two inviscid, non-resistive, 
interacting, magnetized fluids, one charged and 
the other neutral. Since sources and sinks are ex- 
pected to dominate advection in the ion continu- 
ity equation, we treat directly only the continuity 
equation for neutrals, and parameterize the ion 
behavior through equation (2-10). The governing 
MHD equations are thus the equation of continu- 
ity for neutral particles, 

^ + V • (p„v) = 0, (2-la) 
equations of motion for the two species, 

Pn^^ + "^Pn + P„V$G (2-lb) 
+Pnt'm(v„ - Vj) = 0, 

Pi^+VPi+AV$G (2-lc) 

(V X B) X B 

-Pit'm(v„ - Vi) = , 

the induction equation 
r)B 

— = V X (v, X B) , (2-ld) 
and Poisson's equation 

V^^G = 47rGp. (2-le) 

Subscripts i and n denote ions and neutral 
particles respectively. Da/Dat is the convective 
derivative for species a. The ion-neutral colli- 
sion frequency is Vin = Pn{crv)/{mi + m„), and 
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Vni is the neutral- ion collision frequency {piVin — 
Pnt^ni)- In the context of collision rates only, the 
symbol a represents the cross-section for elastic 
collisions; elsewhere it represents surface density. 
The gravitational potential is and v, p, P, and 
B are the velocity, density, pressure, and magnetic 
field, respectively. Wc work on large scales and at 
low temporal frequencies for which the ions and 
electrons are coupled. 

We assume that the ionization fraction in the 
cloud is low. For dense molecular gas which is 
ionized by cosmic rays and recombines on grains, 
n-i/nn ~ Krin^^"^, where if ~ 1.1 x 10~^ (McKee 
et al. 1993) Departures from and generalizations 
of this ionization law are discussed below, see eq. 
[2-10]. Therefore, p„ « p and v„ « v. If the 
neutral-ion collision time is much less than a dy- 
namical time, the ambipolar drift velocity can 
be written in the standard form (Shu 1983): 



(V X B) X B 



(2-2) 



The flat geometry allows significant simplifica- 
tion of the equations by taking a vertical integral 
in the limit of infinitesimal vertical thickness. For 
example, the magnetic force simplifies to: 



lim 

€^0 



dz 



(V X B) X B 



lim —^[BxX + By-y] 



+e 



Bz^h 
27r ' 



(2-3) 



(2-4) 



In the limit of an infinitesimally thin disc or slab, 
the vertical component of the magnetic field B^ is 
continuous with respect to the plane of the slab 
(the z direction), and the horizontal component 
B?i is antisymmetric (reverses sign) with respect 
to the plane of the slab. 

In addition, we assume the vertical component 
of the velocity is negligible compared to the hor- 
izontal components {v^ <C Vx,Vy). Lovelace & 
Zweibel (1997) found that thin disks are generally 
stable to warping modes, so we expect predomi- 
nantly horizontal motion. 

In the limit of zero gas density and thermal 
pressure outside the slab, the external magnetic 
field relaxes instantaneously to an equilibrium 
state, shown in Z98 to be a current-free or poten- 
tial field state. We therefore assume that the mag- 
netic field at 2; 7^ is a potential field. This allows 



us to calculate only the vertical part of the mag- 
netic field in the disc, rather than all three com- 
ponents in a three-dimensional domain, a tremen- 
dous simplification. Limitations of this approx- 
imation, and corrections to it, are discussed in 
more detail in §4.1 and the Appendix. 

The resulting system of 2-dimcnsional equa- 
tions are as follows: The equation of continuity 



da_ 

'di 



+ Vh • {(jwh) = 0, 



of motion 



(2-5a) 



0, (2-5c) 



the definition of the ambipolar drift velocity 

Bz^h 



^Dh = 

the induction equation 
dBz 



2'KVinCri 



dt 



• [{Vh+'^Dh)Bz 



the potential field equations 

„ _ a$A/ 

az 

and Poisson's equation 



dz 



(2-5d) 



(2-5e) 



(2-5f) 
(2-5g) 



(2-5h) 



The gravitational and magnetic potentials are $g 
and $M respectively, a is the surface density, and 
h denotes horizontal components (e.g., voh is the 
horizontal drift velocity). We assume an isother- 
mal equation of state P = cP'u. 

Equations for an axisymmctric disk of small but 
finite half thickness Z were derived by Ciolek & 
Mouschovias (1993). Their equations contain cor- 
rection terms of order Z jP^ where P is the dis- 
tance from the axis of symmetry; these terms in- 
clude magnetic pressure, which provides a restor- 
ing force which we neglect, and corrections to the 
normal direction, which is tilted slightly from the 
vertical because Z depends on P. These terms 
go smoothly to zero in the limit ZjP 0. For 
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canonical values of physical parameters in dense 
molecular clouds, we find Z/R < 1/10. Strictly 
speaking, we should retain the magnetic pressure 
gradient, as Ciolek & Mouschovias (1993) do, be- 
cause it is of the same order as, although gener- 
ally less than, the thermal pressure gradient, and 
also provides a restoring force. However, magnetic 
tension clearly dominates magnetic pressure at the 
long wavelengths of greatest interest here, while as 
the wavenumber increases the rate of ambipolar 
drift increases as well, so that the magnetic pres- 
sure force at short wavelengths is less than it would 
be if the magnetic field were frozen in. Moreover, 
we know that there is a 3D version of the instabil- 
ity which is driven by magnetic pressure alone. We 
study the instability driven by tension, but the ex- 
istence and nature of the instability should be the 
same whether it is driven by tension or pressure. 
For all these reasons, we think that the neglect of 
magnetic pressure is not a major source of error. 

2.2. Nondimensionalization 

All quantities in the problem are scaled by 
a self-consistent set of characteristic quantities. 
Given an initial vertical magnetic field B^o, we 
choose as the characteristic surface density that 
which is marginally stable to collapse in the zero 
temperature limit (Nakano & Nakamura 1978), 



density and magnetic field 



27rGV2 • 



(2-6) 



In a 3-dimensional model of a cold magnetized 
molecular cloud, one logical choice would be to 
use the Alfven velocity as a characteristic veloc- 
ity. The 2-dimensional geometry precludes this 
approach, because the quantity Bzl\/2'Ka which 
arises naturally has dimensions not of (a 2-D 
Alfven) velocity but rather of length^/^/tinie. A 
length scale is thus required. A logical choice 
in this geometry is the scale height of the slab, 
H = a? /2t:(jG, but this is undesirable because the 
problem of most interest is a cold cloud, in which 
the isothermal sound speed — > 0. Instead we 
choose a characteristic length scale L, which will 
be the horizontal domain size, or equivalently the 
largest spatial wavelength in the simulation (Of 
course, L scales out of all final results when ex- 
pressed in dimensional units). The characteristic 
velocity is the Alfven speed for the critical surface 



znaco 

and the characteristic time is simply 

L 

tcO = . 

VaO 



(2-7) 



(2-8) 



The nondimensionalized variables are as follows: 



T 



3 = ^ 



^$G, and 



_ x,y,z 
L ' 



(2-9) 



2.3. Parameterization of Ambipolar Drift 

We assume that the product of the ion surface 
density (Tj and the ion-neutral collision frequency 
is related to the surface density cr ~ cr„ ac- 
cording to the simple ansatz: 



5{(TiVin) 5a 
= a — . 



(2-10) 



If the ionization fraction x scales as the neutral 
density n~', and the scale height H scales as 
(T~^, then a = — 2q. Often q is taken to be 
0.5 (McKee et al. 1993), but a detailed treatment 
of grain dynamics in contracting cores (Ciolek & 
Mouschovias 1994, 1995, 1998) shows that the pa- 
rameter q continuously decreases throughout con- 
traction, and may range from ~ 0.6 to less than 
0.1 as the central density increases by about 6 or- 
ders of magnitude. We have tested the sensitivity 
of our results to the value of a by comparing mod- 
els with a ranging from to 3, and find that the 
results differ by less than 5% as the density per- 
turbations grow from .01 to 1.5 times the mean 
density (This is consistent with linear perturba- 
tion theory, which predicts that the value of a en- 
ters only if the unperturbed initial state has an 
inclined magnetic field; Z98). In view of the in- 
scnsitivity of the results to a. as well as the fact 
that the volume density increases by less than 2 
orders of magnitude in our calculations, we regard 
the ansatz equation (2-10) as adequate. 

To compare with the linear theory, we use a 
nondimensional form of the drift frequency V = 
tc(ikBlQ/2-KCFiVin. (Z98 uses the dimensional T = 
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kBlf^/2'K<7iVin ) In the simulation, it is convenient 
to compute the evolution of the drift indepen- 
dent of spatial wavenumber k. We use the quan- 
tity T/kL, which is r/27r for the lowest spatial 
wavenumber k = 2'k/L. 

The full set of governing equations in conserva- 
tive form are as follows: 



(2-lla) 



d_ 



UJ 



,1 



dr 
kL 



d4>M 



(2-llb) 
(2-llc) 
(2- lid) 
(2-lle) 



Vh-mi^+iyo)], (2-iif) 

r 

kL' 
kL a du) 



(2-llg) 
(2-llh) 



Note that we interpret equation (2-10) as an 
Eulerian relation in equation (2-llh). 

2.4. Physically Reasonable Parameter Regime 

The input parameters for the model are the 
sound speed a, the strength of the initial magnetic 
field relative to the surface density BzQ/2nG^^^(Ta = 
l/ojQ, the drift parameter F, and the constant a 
which determines the perturbation to the collision 
rate (see eq. [2-10]). 

Typical magnetic fields in dense clouds are 
30/xG (Myers & Goodman 1988b), and we choose 
the nondimensionalization length scale (see §2.2) 
to be Ipc, a typical size for a dense cloud or cloud 
core and its close neighborhood. A typical cold 
cloud temperature is lOK, and the average molec- 
ular weight is that of molecular hydrogen with 10% 
helium, m„ = S.QxlO'^^g. This yields the follow- 
ing expressions for the Alfven and sound speeds: 



27rcrcO 



(2-12a) 



B^oLG'/^ 



a 



= 2.4 X 10^° 



ksT 

m, 



~ 0.02 



cm 



lOKj \ LJ \ 
3.9 X 10-24g 



(2-12b) 



(2-12c) 



From now on, will be given in units of v^q. 

The initial surface density is chosen to be ~ 
0.02 g cm"^, or about 100 M0 pe^^. This corre- 
sponds not only to a typical column density for a 
dense cold cloud {Nh ~ 10^^ cm~^), but also to 
the surface density that would result if a spher- 
ical cloud of typical number density (10"* cm"^) 
and typical size (several lO^^cm) collapsed along 
a large-scale magnetic field into a thin pancake or 
disc. The following expression for the normalized 
surface density results: 



Wo 



27rGVVo 



B 



zO 



(2-13) 



1.0 



B,o J VO-02 g cm- 



Thermal pressure raises the critical siirface density 
for gravitational collapse, and for a cloud at T 

' 1 1 \ the value of ujq for the fiducial parameters 
in equation (2-13) is 0.864 of that critical surface 
density Wcrit- 

Thermal pressure imparts a finite scale height 
to the disk 



H = 



2TTaoG 



4 X lO^'^cm 



(2-14) 
0.02 g cm-2^ 



T 
TOK 
3.9 X IQ-^^g 



Clearly, H <^ L. An inclined magnetic field 
exerts a pinching force, compressing the disk and 
reducing H further (Wardle & Konigl 1993). 

Determination of the drift parameter F requires 

an expression for the neutral-ion collision fre- 
quency, ~ 2xl0~^ni cm"^ s"^ (Draine, Roberge, 
& Dalgarno 1983) and the density of ions. We as- 



sume Hi 



1/2 



where K 1.1 x 10"^ 



n. 



1/2 



6 



cm"^/^ (McKee et al. 1993). Using these relations 



(2-15) 



5 X 10* cm 



rir, 



-3\ 2 



27r/L 



The fiducial value of n„ which appears here is 
consistent with the other parameters: n„ = 

Although equation (2-15) shows that F does not 
depend on the scale height H, since the tempera- 
ture of molecular clouds is quite well determined it 
is useful to rewrite F in a way which does depend 
on H and suppresses the dependence on some of 
the other parameters. We have 



F = 



47r-V2(m„G)-''/2 



1/2 




(2-16) 



where we have used the standard values for all the 
constants. Equations (2-14) and (2-16) suggest 
F < 0.1 for typical parameters. 

2.5. Linear Theory 

The collapse rate for a linear perturbation with 
positive thermal pressure can be easily calculated 
(see Z98 for the calculation at T=OK). The physi- 
cal quantities w, Pz, v and governing equations (2- 
5a-2-5h) are linearized and reduced to one horizon- 
tal spatial dimension. Assuming a single Fourier 
mode 



l + ZJ^e^^+^'^^-^l^l, 



(2-17a) 
(2-17b) 
(2-17c) 



leads to the dispersion relation 



7l 



A-1 



■7f 



7-7Sr+7|r = 0, 

(2-18) 



where 7g = ^/cuokL is the nondimensionalized 
gravitational frequency, and jt = akL is the 
nondimensionalized thermal frequency. 

Equation (2-18) has two limits which provide 
some insight into what follows. In the limit F = 
we recover the dispersion relation for waves driven 
by magnetic tension, thermal pressure, and self- 
gravity; the first two forces are stabilizing and the 
last is destabilizing. The system is stable for all 
wavenumbers if wq < 1, but if i^'o > 1, the system 
is unstable for wavenumbers k < fcc(F = 0) = 
H~^{ujq — l)/t^o- (I'^ ^ system of finite size L, k 
is bounded from below by 27r/i, leading to the 
absolute stability criterion wq < i^crit which we 
present below). The maximum growth rate, in 
dimensional form (recall that in equation (2-18), 
7 is given in units of t~Q ), is Jmax(r = 0) = 
nGao{LjQ — l)/(awQ), and occurs at a wavenumber 
fcm(r = 0) = ^kc{T = 0) (in these dimensional 
expressions, a is the dimensional sound speed). 

In the limit of large F, the magnetic field is un- 
coupled from the gas, and the dispersion relation 

reverts to that of an unmagnetized slab. The sys- 
tem is unstable for 7^ > 7|., or A: < fcc(r ^ oo) = 
2wG(Jo/a^ = H~^. The maximum growth rate, 

»), is 



which occurs at km{T 



= ifc,(F 



lmax{y 00) = 7rGao/a. (Again, in these dimen- 
sional expressions, a is dimensional). 

The maximum growth rate, and the wavenum- 
ber at which it occurs, is always less for a magne- 
tized but supercritical cloud than for an unmag- 
netized cloud, and, as expected, the supercriti- 
cal case approaches the unmagnetized case as the 
magnetic fieldstrength decreases to zero. 

In this paper we are interested in clouds which 
are magnetically subcritical, so that they would 
be stable in the limit F = 0, but would be unsta- 
ble if the magnetic field were removed. That is, 
we are interested in clouds with a length much 
larger than the unmagnetized Jeans length. A 
small but nonzero F destabilizes a cloud to per- 
turbations which are stabilized by magnetic fields 
in the absence of ambipolar drift, but would be 
unstable to the Jeans instability in the absence of 
magnetic fields. 

For small soimd speeds, the dispersion relation 
shows the same behavior seen for zero temperature 
in Z98: at low values of the drift parameter F, the 
growth rate of the perturbation 7 is proportional 
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to r. At higher F, however, 7 cx T^^^. As surface 
densities w and wavenumbers k depart from the 
critical values for stability, more ambipolar drift 
is required for the system to show 7 oc T^^^ be- 
havior. To quantify these statements with an ex- 
ample, at Wo = 1-1, = 0.02, F = 0, the critical 
k below which ideal perturbations are unstable is 
k = 9.5455. Very close to criticality, k = 9.6, the 
7 a F^/^ scaling holds for F as small as .001. At 
fc = 10, 7 increases with F faster than T^^^, but 
much slower than linearly, for .001 < F < .01, but 
approaches the V^/^ scaling for .01 < F < .1. As 
F is increased from .001 to .1, 7 increases from 
.378 to 1.75, which is F'^^ scaling. At fc = 47r, 
this scaling has broken down noticeably: as F is 
increased from .001 to .1, 7 increases from .170 to 
1.93, which is F-^'^ scaling. Most of the deviation 
occurs for small values of F; for F between .01 and 
.1, 7 cx F'^^. Thus, for reasonable values of and 
F, the F^/^ scaling law holds quite well even when 
k is as much as 30% below the critical value. 

The addition of thermal pressure increases the 
stability of the disk; more drift (larger F) is re- 
quired for collapse, and more is required to reach 
the transition from 7 cx F to 7 cx T^^^ . An approx- 
imate value for the critical surface density for col- 
lapse, with positive thermal pressure, is obtained 
from solving the dispersion relation (see eq. [2-18]) 
for F=0: 

— = -Tra^ + + TT^a* (2-19a) 

Wcrit 

~ 1 — 7ra^, 
Wcrit ^ l-|-7ro^ (2-19b) 

where the approximations hold for small sound 
speed. Solution of the dispersion relation also 
shows that for a given sound speed and drift pa- 
rameter, there is a single mode with a maximal 
growth rate, and that the modes above a certain 
wavenumber are acoustically suppressed. Figure 1 
shows the growth rate as a function of wavenumber 
for F = 0.1, and most unstable wavenumber, over 
a range of sound speeds. Figure 1 also shows that 
the stability boundary is very near the thermal 
Jeans stability boundary, while the fastest grow- 
ing mode has a much longer wavelength than the 
Jeans wavelength. 

As we will see later, the wavenumber at which 
the growth rate is maximized dominates the struc- 
ture of clumps even into the nonlinear regime. Nu- 
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Fig. 1. — Damping of short- wavelength modes due 
to thermal pressure. The first graph shows the 
linear growth rate 7 of the fastest growing mode 
as a function of spatial frequency / = k/2iT, for 
different sound speeds a^. The second graph shows 
the spatial frequency of maximal growth rate (the 
peak of each curve in the first graph), and the 
highest undamped spatial frequency (where each 
curve in the first graph crosses the x-axis). For 
both figures, F = 0.1, and uJcrit is approximately 
constant: ljq = 0.864(l-|-7ra^). 
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merical solution of the dispersion relation equa- 
tion (2-18) shows that km is always less than 
fcm(r oo), the fastest growing wavenumber for 
the Jeans instability, but also scales with temper- 
ature as 1/T. Therefore, we expect the fragment 
mass to be larger than the thermal Jeans mass, but 
to have the same temperature scaling. This is 
borne out by Figures 1 and 6. 

3. Numerical Simulation 

In order to follow the instability into the non- 
linear regime we have carried out a numerical sim- 
ulation. We use a Fourier collocation pseudo- 
spectral method (Canuto et al. 1988) to solve the 
governing equations. Values of physical quanti- 
ties are stored at discrete points in physical space 
(known as collocation points), and spatial deriva- 
tives are evaluated in Fourier spectral space (hence 
the name "Fourier pseudo-spectral" ) . This partic- 
ular method is well adapted to this problem for 
several reasons. Calculating a spatial derivative 
in Fourier space simply requires multiplication of 
each Fourier component by i times its wavenum- 
ber. All terms involving the horizontal magnetic 
field B?i or gravitational potential are triv- 
ial to evaluate in Fourier space due to the simple 
form of the magnetic and gravitational potentials. 
Nonlinear terms, on the other hand, are trivial 
to evaluate in physical space by simple multipli- 
cation. Finally, growth of different Fourier modes 
can be monitored and controlled explicitly, simpli- 
fying comparison between the nonlinear numerical 
model and the single-wavenumber linear analytic 
results. 

We use a Bulirsh-Stoer time-stepping routine 
with Richardson extrapolation (Press et al. 1986). 
The routine performs several modified midpoint 
method integrations at sub-intervals of the de- 
sired time-step. It then attempts to extrapolate 
to an infinite number of sub-intervals. The rou- 
tine varies the number of explicit (calculated) sub- 
intervals based on the estimated error. In general, 
the full set of governing equations for this problem 
can be integrated with ~ 5 explicit subintervals for 
each time-step oi 5t = 0.1. 

We tested convergence by running the code at 
increasing spatial resolution with the same initial 
conditions. Figure 2 shows the amplitude of a den- 
sity perturbation, computed at different resolu- 



tions, as a function of time. The initial conditions 
had a single wavenumber density perturbation in 
each direction, forming a "checkerboard" pattern. 
(Collapse is described in detail below.) A 16^ grid 
is sufficient to resolve the collapse of such a single 
wavenumber, from the linear regime (exponential 
growth) into the nonlinear regime. Use of a 32^ 
or 64^ grid changes the solution by less than 0.1% 
over most of the time period plotted. Finer grids 
are required to resolve collapse of smaller struc- 
tures, and in runs which contain a spectrum of 
wavenumbers, a 64^ grid was used. 
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Fig. 2. — Convergence of model with increasing 
numerical resolution. Magnitude of the density 
perturbation for the same initial conditions, using 
grid sizes of 4^, 8^, 16^, 322, r^j^jg partic- 

ular rmi had F = 0.1, = 0.1, wq = 0.864 lOcrit 
= 1.135). The oscillations at early times are due 
to the presence of frictionally damped oscillatory 
modes which were present in the initial conditions. 

We have also verified that the code reproduces 
the results of linear theory. If the initial condi- 
tion corresponds to an cigcnfunction of a growing 
mode calculated according to linear perturbation 
theory, with an amplitude of a few percent or less, 
then the disturbance initially grows at the expo- 
nential rate predicted by the linear theory. This 
is shown in Figure 3, which compares the growth 
rates measured from the code (discrete symbols) 
with the continuous curve obtained from solving 
the dispersion relation equation (2-18). The agree- 
ment is generally excellent. 

The model is numerically stable until such time 
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Fig. 3. — Comparison of linear analytical growth 
rate (solid line) with the initial growth rate in our 
simulations (7 when the density perturbation 8lo 
is <10% of mean density w), for a representative 
sound speed = 0.033. Left, the growth rate in 
the simulation (points with error bars) agrees with 
the linear theory (solid line) for a range of initial 
surface density ojo (F = 0.1). Right, the growth 
rate in the simulation is seen to agree with the 
linear theory over a 4 order of magnitude range in 
7. (r varies from 0.01 to 0.3 and wq from 0.3 to 
0.95) 



as power in the higher order wavenumbers grows 
to overwhelm power in the Fourier modes of in- 
terest. At that time, the simulation develops a 
"sawtooth" instability, with large variation be- 
tween alternating collocation points. Power grows 
in these short-wavelength modes from numerical 
noise, whose magnitude is about 10"*^ compared to 
the power in the principal mode (measured in sim- 
ulations whose initial conditions contained a single 
mode). Higher wavenumber modes are also driven 
by the nonlinearity of the problem, and this is the 
dominant physical source of power in those modes. 

Growth of high-order modes can be controlled 
in several ways. Thermal pressure will stabilize 
high order modes, as was seen in Figure 1. In most 
physically reasonable finite cloud temper- 
ature of < 10 K will stabilize the simulation long 
enough to follow the collapse well into the nonlin- 
ear regime. The problem of high-order mode sta- 
bilization is nearly independent of the drift param- 
eter r because an increase (decrease) in ambipolar 
diffusion increases (decreases) the growth rate of 
all modes similarly. Thus the entire physically in- 
teresting part of parameter space is numerically 
accessible and numerically stable in this model. 

3.1. Collapse Rate 

Many runs were computed with initial condi- 
tions corresponding to the eigenfunction of the 
fastest-growing solution of the 3 modes present, 
at each wavenumber, in linear theory. This ini- 
tial perturbation has a single initial wavenumber 
in each direction. Growth initially proceeds expo- 
nentially, with faster nonlinear collapse occurring 
as higher order spatial modes are driven. Non- 
linear behavior is typically seen when the density 
enhancement associated with a perturbation has 
grown to 75%-100% of the mean density. We were 
able to follow the evolution of the system to peak 
surface densities about 10 times larger than the 
mean density (corresponding to a peak volume 
density about 100 times larger than the mean). 

We made a detailed comparison with linear the- 
ory and with simple drift and collapse models by 
studying the early collapse of a perturbation with 
a single spatial wavenumber. Figure 4 shows how 
the growth of a fully nonlinear density perturba- 
tion depends on the drift parameter F. The lower 
three curves describe collapse in subcritical clouds, 
which would be stabilized by the magnetic field in 
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Fig. 4. — Collapse rate for single-wavenumber ini- 
tial conditions, = 0.033, and = 0.954 = 
0.864 Wcrit (lower three curves). These initial con- 
ditions would be stable in the absence of ambipolar 
drift. The solid line is the linear growth rate 7 as 
a function of the strength of the ambipolar drift T. 
The stars connected by a dotted line is the collapse 
rate for the nonlinear simulation. The dashed line 
is the ambipolar drift rate. The pluses connected 
by a dotted line are the collapse rates for the non- 
linear simulation with a high initial surface density 
Wo = 1.278 = 1/0.864 Wcrit- For such supercritical 
collapse, the magnetic field and strength of am- 
bipolar drift should have minimal efi^ects. 



the absence of ambipolar drift. The initial surface 
density (normalized to the vertical magnetic field) 
is ~ 86% of the critical surface density for collapse. 
The growth rate of perturbations in the nonlin- 
ear simulation (calculated from the time for the 
central density of the perturbation to grow from 
1% to 100% of the mean density) is compared to 
the predicted growth rate 7 for linear perturba- 
tions, and to the ambipolar drift rate T. Clearly, 
for the physically expected value of T (^ 0.05 - 
0.10), §2.4), the collapse due to this instability is 
several times faster than simple collapse due to 
loss of magnetic support on a diffusive ambipo- 
lar drift time-scale. For example, when V = 0.1, 
the ambipolar drift rate is 0.1 and the unmagne- 
tizcd collapse rate is 10. The rate of contraction 
found from the simulation is 0.4. At larger values 
of r, the drift rate comes closer to the unmagne- 
tizcd collapse rate. The collapse rate in the simu- 
lation approaches the unmagnetized collapse rate 
because the coupling between the magnetic field 
and the gas is weak. When the drift becomes very 
important (F ~ 10), the drift timescale is very 
short, and the intermediate instability described 
in this paper is less significant. Even for the mod- 
erate sized density perturbations used to create 
Figure 4, (Su) w), the growth rate is larger in the 
nonlinear simulation than in the linear problem, 
especially at relatively large values of F, showing 
that the collapse is accelerated by the nonlinear- 
ity. It is important to note, however, that the 
growth rate shows the same dependence on am- 
bipolar drift in both the linear theory and the non- 
linear simulation: at low F, 7 oc F, and at higher 
F, 70CFI/3. 

The top curve in Figure 4 describes collapse in 
supercritical clouds, in which the magnetic field 
would be insufficient to prevent collapse even if it 
were frozen to the matter. The collapse rate de- 
pends only weakly on the strength of the ambipo- 
lar drift, as expected since the magnetic field is 
dynamically less important. When the ambipolar 
drift strength F becomes large, the drift timescale 
becomes comparable to the collapse timescale, and 
the subcritical and supercritical cases converge. 
Rapid collapse occurs in supercritical clouds due 
to the dynamical weakness of the field, and in sub- 
critical clouds rapid diffusion removes magnetic 
support, quickly rendering them supercritical. 

The evolution of self gravitating, subcritical 
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disks with ambipolar drift was studied previously 
by Ciolek & Mouschovias (1994) and Basu & 
Mouschovias (1995). Ciolek & Mouschovias (1994) 
began with a subcritical (wq = .256), centrally 
condensed equilibrium state - the central surface 
density is 16 times the mean density. Thus, this 
model is more centrally condensed even initially 
than our models are when we terminate the simu- 
lation. The initial ambipolar drift time is 10 times 
the initial free fall time, which corresponds on Fig- 
ure 4 to r ~ 1. The evolutionary timescale in the 
subcritical, quasistatic phase is well estimated by 
the initial ambipolar drift time; after the cloud 
becomes supercritical its collapse rate approaches 
the freefall rate. Although we can extrapolate our 
results to this model only with caution, because 
the initial conditions are so different from ours, it 
does not surprise us that such a subcritical disk 
shows no evidence for the blending of dynamical 
and drift effects that we observe closer to critical- 
ity. 

Basu & Mouschovias (1995) carried out a pa- 
rameter study to determine the effects of the de- 
gree of criticality on the rate of evolution to a 
critical state. They also began with centrally con- 
densed equilibrium models, forming a sequence in 
which the criticality parameter varied from 0.1 to 
0.5. They found that the timescale for evolution 
to the critical state decreased by about a factor 
of 1.5 along this sequence, from somewhat longer 
than the estimated drift time to about 25% shorter 
(another, marginally critical model, collapsed at 
once). Although again a quantitative comparison 
of our models with theirs is difficult because of the 
different initial conditions, it is possible that the 
intermediate contraction rates which they see are 
a manifestation of the coupling between dynami- 
cal and ambipolar drift effects seen in our models. 
It may also be due to the increased central concen- 
tration of the initial equilibrium states along their 
sequence of models. 

The nature and rate of collapse is observable 
in molecular clouds (Evans 1999; Myers, Evans, 
& Ohashi 1999), and an instability with an inter- 
mediate growth rate such as this one can help to 
explain observations that do not fit either of the 
classical scenarios - dynamical collapse or slow dif- 
fusive contraction. Our simulated cloud cores col- 
lapse with slower velocities and on larger physical 
scales than the dynamical inside-out collapse pre- 



dicted when the magnetic field is unimportant, as 
in Shu (1977). Tafalla et al. (1998) and Gregersen 
(1998) have observed cores that appear to have 
such behavior; they find that the regions of inflow 
are too large to fit dynamical inside-out collapse 
models, but that the inflow velocities are too large 
for quasistatic diffusion models. 

3.2. B-p relation 

The correlation between fleldstrength and den- 
sity is an observable quantity which can provide 
insight into the manner in which the magnetic field 
evolves. It is useful to parameterize this relation 
as a power law 

|B|cxp^ (3-1) 

Observations find k ~ 0.5 (Troland & Heiles 1986; 
Crutcher 1999) over several orders of magnitude 
in density. In the case of a highly flattened cloud 
such as we simulate here, it is more convenient to 
define a magnetic field - surface density relation 

|B|occr^. (3-2) 

If the slab scale height H remains constant 
throughout the evolution, then \ = k. If the 
scale height is determined by a balance between 
self-gravity and thermal pressure alone, then an 
isothermal slab obeys H ~ l/tr, and a pH 
implies a ~ p^/^, or k = 0.5A (Crutcher 1999; 
Spitzer 1942). 

If the magnetic field is frozen to the matter 
but not dynamically important, so contraction is 
isotropic, conservation of flux '^mag oc i^|B| and 
mass M oc L^p requires |B| oc p^/^ (k = 2/3). If 
the field is so strong that matter moves one di- 
mensionally, parallel to the fieldlines, then k — > 0. 
Calculations in which a cloud condenses to mag- 
netohydrostatic equilibrium from a uniform ini- 
tial state, with frozen in magnetic fiux of a mag- 
nitude appropriate to the ISM, show anisotropic 
contraction, and the central values of B and p in 
the initial and final states are related by k ~ 0.5 
(Mouschovias & Spitzer 1976; Tomisaka, Ikeuchi, 
& Nakamura 1988). If a cloud is already flat- 
tened and shrinks transversely, M cx L^a, and fltix 
freezing impHes |B| cx (A=l, k=0.5). However, 
rather different input physics leads to a similar 
exponent: simulations of supersonic magnetized 
turbulence, without self-gravity, produce k ^ 0.4 
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if the field is not too strong (Padoan & Norlund 
1999). 

Ambipolar drift generally reduces k below the 
value it would have if the field were frozen in. 
In the models of Fiedler & Mouschovias (1993), 
K averages about 0.2 during the so-called qua- 
sistatic phase. After the quasistatic phase ends, 
the mean value of k is 0.3 as p increases by more 
than 5 orders of magnitude. In the highly flat- 
tened models of Ciolek & Mouschovias (1994), k 
increases smoothly as the cloud evolves from sub- 
critical to supercritical, reaching peak values be- 
tween 0.4 and 0.5 and being about 0.3 at critical- 
ity. In our simulations the magnetic field - sur- 
face density exponent A varies between 0.35 and 
0.65. As shown in Figure 5, the exponent de- 
creases as the central density of a clump increases. 
As collapse proceeds, not only does the density 
increase, but the magnetic field curvature also in- 
creases as fieldlines are dragged into the conden- 
sation. Both effects increase the ambipolar drift 
velocity = Bz^h/inavni and thus the rate 
of flux loss from the clump. Our models do not 
show the increase of A toward its frozen flux value 
as the central density increases seen in Ciolek &; 
Mouschovias (1994), because we follow only the 
early stages of contraction, in which the velocity 
is well below the frccfall value. The exponent A 
also decreases as the amount of ambipolar drift F 
increases, as would be expected, and as the cloud 
temperature increases. The latter effect results 
from the decreased efficiency of this instability in 
warm clouds. The collapse rate 7 decreases with 
increasing temperature as shown in §2.5, and the 
collapse time is longer relative to the ambipolar 
drift time, so more flux can leak from the clump 
as it collapses. 

Comparison with observations depends on the 
assumed disc scale height H and the central densi- 
ties of observed cloud cores. Crutcher (1999) finds 
K = 0.47 for cloud cores with densities lO^'^cm"*^ 
< uh ^ 10 ''■^cm"'^. His data are strongly weighted 
by observations which only measure upper lim- 
its for the magnetic fieldstrength, and omission 
of those data points results in an exponent of k 
= 0.3 over the same range of densities. We mea- 
sure 0.3 < A < 0.7 in simulated cores. If the slab 
scale height H is constant, then 0.3 < k < 0.7, but 
if the slab obeys H ~ 1/a, then 0.2 < k < 0.35. 
Our simulations are thus consistent with the small 
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Fig. 5. — The magnetic fleld - density exponent 
A. (|B| cx (T^). The exponent decreases (increased 
flux loss from a collapsing clump) with increas- 
ing central density Sui, increasing ambipolar drift 
strength F, or increasing cloud temperature T. 
The first graph shows variation with 6u), the sec- 
ond variation with F at constant = 0.05, and 
the third, variation with T at constant F = 0.1. 
In all simulations, loq = 0.864 iVcrit- 
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number of available observations. 

3.3. Clumps and Fragments 

Real molecular clouds have density variations 

on many scales, or a spectrum of spatial wavenum- 
bers. The linear theory for this instability (§2.5) 
predicts that a single mode will dominate collapse. 
Fragments of a single mass will form, with that 
mass depending on the temperature and degree of 
criticality, with little dependence on the strength 
of ambipolar drift, provided that it is small. We 
have simulated many (~ 100) clouds in which the 
initial density perturbations have a broad spatial 
Fourier spectrum and random phase. The real 
part of the initial Fourier spectrum of the den- 
sity perturbation is a Gaussian centered on fc = 
(with the omission of the fc = mode itself). 
The FWHM is 0.33 

kmax^ where kmax is the high- 
est spatial mode in the simulation, so a range 
of low-order modes have similar initial strength, 
and there is significant initial power even in some 
acoustically damped (high k) modes. After some 
time, the clouds coalesce into a small number of 
fragments, each of which we define to be a re- 
gion with a > ao- The size of the fragments is 
well-predicted by the linear theory. For example, 
when = 0.02, the fastest linearly growing mode 
is fc/27r ~ 2.5, and in the simulation, the modes 
{kx, ky) /2Tr = (1,3) and (3,1) have much larger am- 
plitudes than other spatial modes. 

Our simulations produce a wide variety of 
clumps and fragments, as expected for systems 
with random initial density fiuctuations. Figure 
6 shows the fragment masses at a fairly early 
stage of collapse (the density perturbation is ^ 
50% greater than the mean density). There is 
considerable scatter, but the linear theory is a 
good guide to the average fragment mass. The 
clump masses are similar to the masses of the 
supercritical cores formed in some previous sim- 
ulations (Fiedler & Mouschovias 1993; Basu & 
Mouschovias 1994; Ciolek & Mouschovias 1994; 
Ciolek & Konigl 1998). 

At sufficiently large T, the spatial mode with 
the largest linear growth rate is the fundamen- 
tal mode in the simulation domain (A = L), and 
fragments larger than this cannot form in a peri- 
odic simulation. This explains the apparent flat- 
tening of the clump mass - temperature relation 
seen at the higher temperatures, as clump masses 



are bounded above by the mass of that "lowest 
mode" . 
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Fig. 6. — Masses of cloud fragments as a function 
of temperature. Masses in the simulation (crosses) 
follow the mass of the fastest growing mode in lin- 
ear theory (solid line) unless that mode exceeds 
the size of the simulation, (see text) These runs 
have r = 0.1, cjo = 0.864 uicrit- Masses from the 
simulation are plotted assuming L = 1 pc, but the 
agreement between the mass predicted by linear 
theory and the mass in the simulation is indepen- 
dent of L. 

The question of whether there is a characteristic 
mass for molecular cloud substructure is an impor- 
tant one. Several observational studies have found 
a power law distribution of clump masses over sev- 
eral orders of magnitude, dN/dM oc M~p, where 
p ~ 1.5 - 1.7 (Blitz 1993). Kramer et al. (1998) 
present evidence that the power law extends far 
below 1 Mq. However, in the Taurus molecular 
cloud there is evidence of a minimum scale of a 
few tenths of a parsec, corresponding to several 
solar masses (Blitz & Williams 1997). Good- 
man et al. (1998) have argued for an inner scale 
of 0.1 pc, which they identify with a transition to 
what they term "velocity coherence" . These inner 
scales are of the same order as the thermal Jeans 
length, and also close to the cutoff wavelength be- 
low which Alfvcn waves are critically damped due 
to strong ambipolar drift (McKee et al. 1993). 
Our results suggest that there is another scale, 
which is somewhat larger, in magnetically subcrit- 
ical clouds with weak ambipolar drift. 
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3.4. Velocity Structure 



Unlike axisymnictric; c;ollapsc rnodcils, thcisc 
simulations allow the study of asymmetric col- 
lapse, clumps with complicated morphology, rela- 
tive motion of clumps, and their internal velocity 
and vorticity fields. We find that collapse is often 
asymmetric and that significant vorticity is gen- 
erated by the instability (although of course the 
net angular momentum remains zero in our simu- 
lations; its absolute value is more than an order of 
magnitude smaller than the estimated numerical 
errors). 

Figure 7 shows four examples of collapsing frag- 
ments. The left panels show contours of constant 
surface density, arrows indicating the local direc- 
tion and magnitude of velocity, and bold arrows in- 
dicating the center of mass motion of each clump. 
The velocity field shows infall towards clumps, but 
there is also a visual impression of "swirling" or 
rotational motion. This is borne out by the right 
panels of Figure 7, which show contours of con- 
stant vertical vorticity overlaid on contours of con- 
stant density. We show below that the magnetic 
field generates local vorticity. 

The velocity field indicates that the collapse is 
in most cases asymmetric, with e.g. much greater 
infall velocities on one side of the clump than the 
other, and significant nonradial motion. Clump 
mergers are possible - one is quite likely taking 
place in the bottom panel of Figure 7 - but the 
bulk motions are significantly slower than the in- 
fall velocities internal to clumps. The reverse is 
generally true in molecular clouds {e.g. Blitz 
1993), and this result in our simulations is a man- 
ifestation of the fact that the velocity field in the 
system is weak. We return to this point in §3.5. 
Sometimes the clumps move apart, but this is am- 
biguous in a periodic domain. The infall veloci- 
ties within individual clumps are of order 0.25 a, 
and are more consistent with observations (Evans 
1999; Myers, Evans, & Ohashi 1999). 
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Fig. 7. Velocity, density, and vorticity for some typical simulations. 

The left panel of each row shows the velocity field and density contours. Contour levels denser than the 
mean density are solid, those less dense than the mean are dashed. The velocity scale is indicated by the 
arrow in the lower left hand corner, whose length is 25 m/s. The density- weighted mean velocity of each 
clump is indicated with a thick arrow, and the scale is 10 times the general velocity scale (the arrow in the 
lower left would represent 2.5 m/s) 

The right panel shows contours of vertical vorticity (V x v)^ for the same cloud. Positive (solid) and negative 
(dashed) vorticity are plotted along with the surface density (dotted). 
For all runs, = 0.025, F = 0.1, loq = 1 = 0.927 Ucrit- 
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Fig. 7. — continued. 
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Visual inspection of Figure 7 shows that the 
clumps are distinctly noncircular and quite elon- 
gated in shape. Accounting for the third dimen- 
sion, our clumps should be considered triaxial or 
prolate. These shapes are consistent with the mea- 
sured and inferred shapes reported by Myers et 
al. (1991), Rydcn (1996), and Ward-Thompson, 
Motte, & Andre (1999). 

Figure 8 shows a cross section of the surface 
density and velocity profiles across the short axis 
of one particular collapsing clump. The density is 
much more peaked than the velocity at this early 
stage of collapse; the FWHM of the infall speed is 
several times larger than the FWHM of the den- 
sity peak. Both the velocity and density profiles 
are clearly, but not grossly, asyimuctric. Although 
it is premature to compare these density and ve- 
locity profiles with observations of infall, it is en- 
couraging that we see evidence for extended in- 
ward motions as have been reported (Tafalla et al. 
1998; Gregersen 1998; Evans 1999; Myers, Evans, 
& Ohashi 1999). A general feature of infall onto 
a line mass such as a filament or strongly pro- 
late object is that the velocity decays more slowly 
with distance from the mass centroid than for in- 
fall onto a spherically symmetric, centrally concen- 
trated object. This may be the main effect that 
produces the extended infall. 

The vorticity generated in swirling motions is 
V X f ^ 5/tco ~ 8/Myr. Collapse proceeds on 
a timescale of several Myr, so the swirling and 
rotation of clumps is not insignificant, although 
it is dominated by infall and we have not found 
evidence for clumps torn apart by shear. Magnetic 
braking, which is excluded from our calculation by 
the potential field approximation, would reduce 
rotation. We estimate the magnetic braking rate 
in §4.1. 

Figure 7 shows that the vorticity maxima are 
displaced from the density maxima. In order to 
understand this, we derive an evolution equation 
for the z component of vorticity, ojz , by taking the 
curl of the equation of motion (2-5b) 

+ V,. • {u,Vh) = Bft X ^h^- (3-3) 

According to equation (3-3), the generation of 
vorticity is second order in the amplitude of the 
fluctuation. There is generation of vorticity to 
flrst order only if there is a zero-order inclined field 
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Fig. 8. — Cross section of density perturbation (5a; 
(solid) and infall velocity (dotted) across a clump. 
Infall velocity increases towards the center of the 
clump (see text). This run is the third example 
in Figure 7 - F = 0.1, = 0.025, wo = 1- = 
0.927 ojcrii- The collocation points are separated 
by about 5 x 10^^ cm in this run, so the density 
and velocity profiles are quite well resolved, except 
at the very center of the density peak. 

(Z98). We can imderstand the spatial pattern of 
vorticity as follows. Despite ambipolar drift, the 
contours of constant track the contours of con- 
stant (7 quite well. Therefore, the gradient of Bz fa- 
is maximized toward the outer edge of a clump, 
not at its center. Note that if the magnetic flux 
were perfectly frozen to the matter, B^/cr would 
retain its initial constant value and there would 
be no vorticity production at all. However, real 
clouds probably have spatially varying Bz/cr, so 
in general, vorticity production does not require 
ambipolar drift. The maximum of B/j, like the 
maximum gradient of Bz/ (t, is displaced from the 
clump center. Equation (3-3) shows that therefore 
vorticity is generated off-center as well, and gener- 
ally changes sign across the clump, so that clumps 
are associated with vortex pairs. Moreover, equa- 
tion (3-3) shows that an axisymmetric clump does 
not generate vorticity. The dynamical pressure of 
the vortices accentuates the non-axisymmetric na- 
ture of clump contraction, and appears as stream- 
ing motions along the major axis of the clump. 
Although in principle equation (3-3) suggests that 
the vortical velocity is scaled by the Alfven speed, 
in our numerical models the vortical velocities are 
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rather small, somewhat less than the infall veloc- 
ities. This is large enough to noticeably elongate 
the clumps, but not enough to tear them apart by 
shear. 

3.5. Energy Redistribution 

Z98 suggested that this magneto-gravitational 
instability might generate significant turbulent ki- 
netic energy by releasing energy contained in the 
background magnetic field. Analysis of the total 
gravitational, magnetic, and kinetic energy (Fig. 
9) in these simulations shows that the absolute 
values of all three forms of energy grow exponen- 
tially during collapse. The magnetic energy, which 
is the fluctuation energy integrated over the space 
outside the disc, dominates at all stages of col- 
lapse in these clouds, which would be stabilized 
by the magnetic field in the absence of ambipo- 
lar drift. The initial magnetic field is uniform in 
these simulations, and so there is no stored mag- 
netic energy available for conversion to turbulent 
motions. The relatively low kinetic energy in the 
models is a signature of the importance of diffu- 
sive, as opposed to dynamical, effects. In ideal 
MHD turbulence with self gravity one would ex- 
pect cquipartition between the kinetic and poten- 
tial energies (Zweibel & McKee 1995). In this case 
the potential energy is the sum of the gravitational 
and magnetic energies, but Figure 9 shows that the 
kinetic energy is about an order of magnitude less 
than the equipartition value. 

4. Discussion of Approximations 

It is difficult and computationally expensive to 
simulate the evolution of magnetized, self gravitat- 
ing molecular clouds in three dimensions with suf- 
ficient resolution to capture all the relevant phys- 
ical processes. In this study we focussed on the 
two-dimensional instabilities of a sheetlikc cloud 
surrounded by a conducting medium without in- 
ertia. This allowed us to approximate the external 
magnetic field as current free, and to work in only 
two spatial dimensions. In the next two subsec- 
tions we discuss the accuracy of these approxima- 
tions. 

4.1. Potential Field Approximation 

In the Appendix, we show how to extend each 
Fourier component of the magnetic and velocity 
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Fig. 9. — Total energies in the simulation domain. 
Magnetic (dashed) and kinetic (dotted) energies 
are plotted with the absolute value of the gravi- 
tational (solid) energy. The absolute values of all 
forms of energy increase exponentially as collapse 
progresses, but the relative distribution does not 
change significantly (see text.) For this particular 
run, r = 0.1, ujo = 0.973 = 0.864 ujcmt, and = 
0.04. 

perturbations above and below the sheet. Al- 
though the perturbations become highly nonlin- 
ear within the cloud, nonlinear effects outside 
the cloud are weak as long as the velocities in- 
side the cloud are sub-Alfvenic with respect to 
the intercloud medium, which is expected for low 
ambient density. The compressive part of the 
cloud velocity field generates evanescent, fast mag- 
netosonic waves which decay exponentially away 
from the disc, while the vortical part generates 
Alfven waves which propagate away from the disc. 

Both the magnetosonic and Alfven waves 
slightly change the horizontal magnetic field per- 
turbation in the disc, thereby changing the mag- 
netic force from its value in the potential approxi- 
mation. If we define an external Alfven timescale 
TAe = {kvAe)~^ for wavcnumber k in the disc, 
and let the timescale for the perturbation in the 
cloud be Tc, then according to equation (AS) the 
correction to the force due to compressive motion 
is of order (tab/tc)^, while the correction due to 
vortical motion is of order (T/ie/T"c)Co/I?o, where 
Cq/Vq is the ratio of the amplitude of compres- 
sive to noncompressive motion. If the motions in 
the disc were Alfvenic, tab/tc would be of order 
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the ratio of the cloud density to external density 
{Pc/ pe)^^"^, but the perturbation frequency is sub- 
Alfvenic, so the ratio of timescales is even larger. 
Moreover, the perturbations arc primarily com- 
pressive rather than vortical. Thus, the error in 
the force incurred by the potential approximation 
is likely to be small. Even if the density contrast 
were only 10^, and the motions in the disk were 
Alfvenic and purely vortical, the potential field 
approximation would still be accurate to 10%. 

The Alfven wave flux tends to suppress the in- 
stability, and removes vorticity from the cloud, at 
a rate that we can quantify. We define an energy 
damping time 7^ as the ratio of outward propagat- 
ing energy flux to the vertically integrated wave 
energy in the disc. Prom equation (A9), 



Id = 



2peVAe 



/72 



+ 



(4-1) 



If we replace a by the critical surface density 
assume that the vertical flelds in- 
side and outside the cloud are the same then 7^ is 
just the gravitational frequency for the intcrcloud 
medium, reduced by the ratio of vortical to total 
kinetic energy 



7d = (47rGpe) 



1/2. 



r2 



Co' I + I T^l 



(4-2) 



Since self gravity is presumably negligible in 
the low density intcrcloud medium, the energy 
loss rate is negligible as well. Loss of vorticity 
is measured by the magnetic braking rate 7^6, 
which can be shown by a similar argument to be 
7^6 - {AtiGp,Y/\ 

We thus see how the potential field limit is 
approached as the density contrast between the 
cloud and intcrcloud medium increases. At the 
late stages of clump formation the potential field 
becomes highly distorted and develops partially 
closed topology, but we can ignore this for the rel- 
atively mild density contrasts studied in this pa- 
per. 

4.2. Approximations to the Gas Physics 

The two dimensional approximation has a long 
and venerable history in galactic dynamics and ac- 
cretion disc theory, as well as in studies of molec- 
ular clouds, and its errors for self gravitating sys- 
tems are reasonably well understood. We expect 



the approximation to be reasonably good as long 
as the clump diameters exceed the disc thickness. 
The instability discussed in this paper has a 3D 
analog (Z98), but it must be treated by other 
means. 

We assumed that the gas has an isothermal 
equation of state. This is a reasonable description 
of the kinetic pressure - density relationship, be- 
cause of the high radiative efliciency of molecular 
gas. However, if the pressure were due to unre- 
solved turbulence the medium would generally be 
less compressible; for example, Alfven wave pres- 
sure follows density according to a 3/2 law (McKee 
& Zweibel 1995). This would make the medium 
more stable by increasing the value of 7t (see eq. 
[2-18]), as would retention of magnetic pressure. 

We took a uniform sheet at rest as an initial 
condition. This has the advantage of simplicity, 
but it means that there is no free energy stored 
in the background magnetic field. Thus, we have 
not tested the conjecture that the instability can 
convert magnetic free energy to turbulent energy, 
which was proposed in Z98. 

We have treated ambipolar drift in the strong 
coupling approximation, and have implicitly as- 
sumed that vo < 20 km s"-'^ (otherwise the rate 
coefficient would change). This is reasonable as 
long as the ion-neutral collision time, which is of 
order 5 x 10^n~^ s, is less than other timescales in 
the problem. 

We have parameterized the relationship be- 
tween the collision rate and the surface density 
by an exponent a. In order to do better we would 
need a three dimensional model of the sheet and 
might need to follow the ionization as well. The 
results of linear theory are rather insensitive to 
the value of a, which suggests that it need not be 
calculated very accurately in the present models. 

5. Summary 

The study of axisymmetric contraction of 
weakly ionized, self gravitating, magnetized clouds 
has proceeded quite far (Basu & Mouschovias 
1994; Ciolek & Mouschovias 1993, 1994; Safier, 
McKee, & Stahler 1997; Ciolek & Konigl 1998). 
In this paper, our emphasis has been on the ini- 
tial breakup of a cold magnetized cloud gas into 
fragments and the early stages of their magnetic 
flux loss and contraction. We include self gravity. 
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magnetic tension, and ambipolar drift, but we do 
not include detailed chemistry of grain physics, 
choosing instead a simple parameterization of the 
ionization. Wc follow the evolution for a shorter 
time than the isolated cloud collapse studies, but 
we impose no symmetry constraints or initial den- 
sity or velocity structure. 

We study highly flattened clouds, with an ini- 
tially perpendicular magnetic field, which are 
slightly subcritical. The linear theory of collapse 
in such geometry (Z98, T=0; this work, T>0) 
predicts collapse on an intermediate timescale, 
faster than the; diifusive timescale set by ambipo- 
lar drift, but slower than the dynamical timescale 
of free-falling inside-out collapse. The linear the- 
ory also predicts the existence of a single spatial 
wavenumber with maximal growth rate, with suf- 
ficiently short wavelengths stabilized by thermal 
pressure. This naturally suppresses power at short 
wavelengths, which is important for the success of 
the spectral method we employ in the simulations. 

We simulate collapse in clouds with random ini- 
tial density perturbations which grow from <0.1% 
of the mean density to 5-10 times the mean den- 
sity. We confirm the intermediate collapse rate 
predicted by linear theory (§3.1), although the 
nonlinear collapse rate is faster than the linear 
rate. These intermediate rates are consistent 
with some recent observations of infall in molecu- 
lar clouds (Evans 1999; Myers, Evans, & Ohashi 
1999). 

Wc show that clouds fragment into clumps with 
size corresponding to the wavelength of the spatial 
mode of maximal linear growth rate (§3.3), gener- 
ally 1-10 Mq. Collapse is asymmetric and com- 
plex (§3.4), and generally forms prolate clumps, 
for which there is observational evidence (Myers 
et al. 1991; Ryden 1996; Ward-Thompson, Motte, 
& Andre 1999). Sometimes the clumps are in mu- 
tual orbit, although the typical clump separation, 
a few tenths of a parscc, is too large to be rele- 
vant to the formation of binary stars. The mag- 
netic field drives the growth of local vorticity, typi- 
cally in the form of vortex pairs which straddle the 
clumps and are associated with streaming motions 
along them. 

Considerable magnetic flux is lost from the col- 
lapsing clumps, consistent with the observation- 
ally determined |B| a (Troland & Heiles 
1986; Crutcher 1999) (see also §3.2). This flux 



loss is consistent with other calculations of cloud 
evolution (Fiedler & Mouschovias 1993; Ciolek & 
Mouschovias 1994) (although ambipolar drift is 
not necessary to bring about this relationship, 
either for isolated (Mouschovias 1976) or turbu- 
lent, highly structured (Padoan & Norlund 1999) 
clouds). As magnetic flux is lost and the surface 
density increases in the central regions of a con- 
tracting core, further fragmentation might ensue. 

One prediction of the linear theory, namely that 
the instability could convert magnetic free energy 
to turbulence, has not been borne out by the simu- 
lations. This may be due to the fact that the initial 
magnetic field is completely uniform and therefore 
carries no free energy. Although significant mag- 
netic curvature develops late in the runs, the cloud 
has already become quite dynamical. This predic- 
tion awaits future tests with a more stressed initial 
state. The relative motions of the clumps shown 
in Figure 7 are about an order of magnitude less 
than the relative motions of clumps separated by 
a few tenths of a parsec in real clouds (Goodman 
et al. 1998), although the infall velocities in the 
simulation are comparable to measured velocities 
(Myers, Evans, & Ohashi 1999). 

An interesting area of future work would be 
to extend this study to true three-dimensional 
clouds. The linear theory (Z98) indicates that this 
instability exists in three as well as two dimen- 
sions, and that the growth rate is still intermedi- 
ate to slow diffusive contraction and fast dynam- 
ical collapse. Construction of a nonlinear model 
in three dimensions would be more difficult than 
the two dimensional models developed here, but 
could prove interesting. In this vein, we find the 
recent successful fit of observations of L1544 with 
a nearly critical model (Ciolek & Basu 2000) en- 
couraging. 
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5-4063 to the University of Colorado, as well as 
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A. Appendix 

In this Appendix wo drop the potential field approximation and calcnlate the; response of the intereloud 
medium to motions within the cloud. This allows us to estimate the errors incurred by assuming a potential 
field. 

We carry out the estimate using linearized, ideal MHD theory. This is more accurate in the intereloud 
medium than it would be in the disc, because the intereloud or external Alfven speed VAe is relatively 
high while ion-neutral friction is weak. We assume the equilibrium intereloud field is uniform and vertical 
(Bq = zBq). We choose to work in the half space z > (the results are similar in the other half space). In 
linear theory, the motions are purely horizontal, and we can derive the following pair of decoupled equations 
for the divergence V and vertical component C of the curl of the velocity 

(^-4eV'j© = 0, (Ala) 

Equations (Ala) and (Alb) represent fast magnetosonic waves and Alfven waves, respectively. In general, 
both types of waves are generated by the motions in the cloud. 

In order to make progress, we assume plane wave horizontal behavior and exponential behavior in time, 
so that all perturbations depend on {x, y, t) as exp(7t + ikxX + ikyy), where 7 may be complex: j = iuj + u 
with both LJ, f > 0. Then 

C = i{kxVy - kyVx); V = {{k^Va, + kyVy) (A2) 

can be calculated at z = in terms of the motions on the disc. The vertical extensions of these quantities can 
be found from equations (Ala) and (Alb), choosing outward going or exponentially decaying wave solutions. 
For the Alfvenic part, 

C = Coe-'>'^'; kA=—, (A3) 

VAe 

where Co is the value of C at 2; = 0. For the magnetosonic part. 



P = Poe-*^"^; fcM^fc^l I + T2V ) ' (A4) 

where k\ = k"^ + ky. In equation (A3) we have written the vertical dependence as a propagating wave, and 
in equation (A4) as an evanescent wave. Although both kA and ku are complex, because 7 is complex, our 
notation reflects the salient aspects of their behavior. The magnetosonic wave is almost purely evanescent 
because the wave frequency is much less than the disc Alfven frequency, which in turn is much less than the 
intereloud Alfven frequency. The Alfven wave has a substantial propagating component and decays in space 
as long as the disturbance is growing in time, which is purely a result of causality. 

We now calculate the perturbed magnetic field components at the disc. According to the linearized 
induction equation, 

5(5B± (9v_L 

^ = -BoV. (A5b) 
Inverting the definitions of C and T) for the velocity components gives 

Vx = -rri^y^ ~ kx'D); Vy = -^{k^C + kyV). (A6) 



2 
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Using equations (A3), (A4), and (A6) in equations (A5a) and (A5b) gives the field components at ^; = 

SB^ = -^{kAkyCo + ikMkxVo), (A7a) 

K_l_7 

SBy = -^{-kAkxCo + ikMkyVo), (A7b) 

«;_i_7 

SB, = -^Vq. (A7c) 
7 

We can use equations (A7a-A7c) to compare the MHD solution with the potential field limit. The ratios 
of the perturbed horizontal to vertical field components at ^; = can be written as 

^ = -ik. ( 1 + jA-] + X ^^)r^^- m 

SB, y k\v\^j k±VAeT>o 

Equation (A8), together with equation (A3), shows that in the limit vab oo the potential field solution is 
exact. 

The Alfvenic part of the disturbance, as a propagating wave, removes both energy and angular momentum 
from the cloud. The energy flux (accounting for both kinetic and electromagnetic energy, and for waves 
propagating in both directions away from the disc) is 

= 2pe^-r§r^VAe- (A9) 

Kj_ 

In §4.1 we use equation (4-1) to derive the rate at which the perturbation in the disc is damped by outgoing 
waves. 
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